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We present a modal approach to calculate finite temperature Casimir interactions between two 
periodically modulated surfaces. The scattering formula is used and the reflection matrices of the 
patterned surfaces are calculated decomposing the electromagnetic field into the natural modes of 
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the modal approach and compared to experiment for validation. The Casimir force from a two 
dimensional periodic structure is computed and deviations from the proximity force approximation 
examined. 
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I. INTRODUCTION 

The dynamics of a classical or a quantum field dras- 
tically depends on the external boundary conditions im- 
posed on it. These boundary conditions lead to a modifi- 
cation of its power spectrum with consequences on mea- 
surable quantities like the radiation pressure. One of the 
most notable examples of this kind of phenomena is the 
Casimir effect [1 j: in its original formulation the change 
in the spectral density of zero-point fluctuations of the 
quantum electromagnetic field induced by the presence of 
two perfectly reflecting parallel plates turns into a net ra- 
diation pressure that pushes one plate towards the other 
or, in other words, into an attractive force between the 
plates. 

Recently, thanks to technological advancements, we 
have witnessed an increased interest in Casimir interac- 
tions [2H9]. Indeed, the Casimir force offers new possibil- 
ities for nanotechnology, such as actuation in micro- and 
nanoelectromechanical systems (MEMS and NEMS) me- 
diated by the quantum vacuum. However, it also presents 
some challenges since the same force is generally recog- 
nized as one of the possible sources of stiction and conse- 
quently of malfunctioning of these devices. Researchers 
have therefore started an intense theoretical [TQHT7] and 
experimental [THJ [19] program in order put some the- 
oretical constraints [T0HT3] and to understand how to 
engineer the strength and possibly also the sign of the 
Casimir force - a repulsive Casimir force would provide 
an anti- "stiction" effect. Inevitably a lot of attention 
has been focused on the role of the boundary conditions 
and very recently on the interplay of material properties, 
temperature, and geometry. 

In this paper we present our results for the computa- 
tion of finite temperature Casimir forces between periodic 
nanostructures using a modal approach. Our calculation 
is based on the scattering approach to the Casimir ef- 
fect which gives the Casimir force between two objects 
starting from their scattering properties. In our case this 
reduces the problem to the calculation of the reflection 
matrices of the periodic nanostructures. We will pay par- 
ticular attention to 2D lamellar gratings, which are peri- 



odic metallic and/or dielectric structures that consist of 
planar layers. Many complex 3D structures, such as pho- 
tonic crystals and metamaterials, can be thought of as 
being constructed from individual 2D extruded metallic 
and dielectric strata [20]. The problem of the reflection 
of an electromagnetic field impinging on a periodic struc- 
ture is a topic that has a huge literature, and there are 
a large variety of methods to efficiently compute the re- 
flected and transmitted electromagnetic power from such 
a surface (for a review see, for example, [2lJ|22]). Among 
the most famous methods one finds the differential ap- 
proach [23] , the integral approach [24] , and the modal ap- 
proach [20]. Our numerical results are based on the last 
technique that, for its characteristic of simplicity, flexi- 
bility, and stability lends itself to be the most adequate 
to our purpose. (See [25] for a review of several numeri- 
cal techniques specifically adapted to the computation of 
Casimir forces.) 

While other more general numerical scattering tech- 
niques exist, the modal expansion of the electromagnetic 
field provides insight into the anatomy of the Casimir 
force. This microscopic analysis provides an understand- 
ing of the nature of the electromagnetic fluctuations that 
give rise to the Casimir force and a means to modify 
their contribution to the Casimir force [26H30] . These 
dominant eigenmode and frequency contributions to the 
Casimir force can be identified within the modal expan- 
sion and their influence under perturbation of the per- 
mittivity of the structure explored. Our modal expan- 
sion, based on a planewave expansion of the fields and 
a Fourier decomposition of the permittivity of the struc- 
ture, is limited to periodic structures but is the natural 
choice to examine the Casimir force in lamellar structures 
like photonic crystals and metamaterials. 

The paper is organized as follows. In Section [IT] we 
briefly review the calculation of the Casimir interaction 
free energy and force at finite temperatures within the 
scattering approach, and we specialize the general for- 
mula to periodic structures. In Section III we describe 



the modal approach and its application to generic (non- 
lamellar) periodic structures. Our finite-temperature nu- 
merical code is then benchmarked in Section IV against 
experimental data for the Casimir interaction between 
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a metallic sphere and a doped Si ID lamellar grating 
[31]. Related zero-temperature computations using the 
differential approach have been performed in [32H34] for 
dielectric ID lamellar gratings, and in [35] for metallic 
ID sinusoidal gratings. We also compute Casimir forces 
in two-dimensional structures consisting of 2D array of 
silicon pillars, or a 2D array of square holes (the com- 
plementary structure to the array of pillars). In the last 
section we discuss our results and prospects for future 
studies. 



II. THE CASIMIR FREE ENERGY 
A. General framework 

In the last few years several authors have devel- 
oped many powerful (semi) analytic [36-43 and numeri- 
cal [44j [45] approaches to calculate Casimir forces. They 
allow for the treatment of quite general geometric and 
material setups going beyond the simple plane-plane ge- 
ometry of Casimir's seminal calculation. For our pur- 
poses, the scattering approach [39] [42! is the most con- 
venient formulation. The advantage of this technique is 
that it exclusively relies on knowledge of the reflection 
properties of the objects seen as isolated scatterers for the 
electromagnetic radiation. Indeed, for linear magneto- 
dielectric media, i.e. media that are completely char- 
acterized by permittivity ^(oj) and permeability ^(cj) 
tensors, virtual and real photons are treated at the same 
level [46 , so that vacuum fluctuations are scattered in 
the same way as real fields. 

The interaction free energy between two bodies in ther- 
mal equilibrium at a temperature T is given by 

-j CO 1 

Ho) =gYl Tr lo S I 1 - • X M ■ ^2 • X 21 (a)] , (1) 

where f3 = l/k&T. Si is the scattering operator charac- 
terizing the i-th object seen from the other object as if it 
were isolated. Xi2(a) and X2i(a) are the so called trans- 
lation operators [39] and they carry no information about 
the objects but only depend on the distance a between 
the origins of appropriately chosen coordinate systems in 
each body. All operators depend on / through the Mat- 
subara (imaginary) frequencies uj\ = i^i = i2irl/hf3 (the 
prime in the summation indicates that the I = term is 
weighted by 1/2), and their detailed expressions depend 
on the functional basis we choose to describe the elec- 
tromagnetic field. The symbol "Tr" indicates the trace 
over all spatial degrees of freedom, making the final result 
basis independent. 

B. Periodic systems 

The expression ([!]) , being the infinite dimensional trace 
of a very involved operator, is in general an extremely 
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FIG. 1: Casimir set-up for two periodic structures (possibly 
mult i- layered) parallel to each other and separated by a gap 
a. 



complicated object. Most of the difficulty consists in 
finding the appropriate functional basis in which the op- 
erators in can be efficiently calculated and, when pos- 
sible, reduce to simple expressions. Most of the time the 
suitable basis for Si is completely different from the one 
for c>2 • In such circumstances other matrices representing 
the change of bases have to be included in the represen- 
tation of ([!]) (see, for example [(39) [47]). 

Fortunately, the symmetry properties of periodic sys- 
tems suggest specific bases that allow for further manip- 
ulation and useful simplifications. This is the case of the 
simple setup where two gratings are facing each other 
shown if Fig. ([I]). The periods of the gratings are as- 
sumed to be the same for simplicity (for gratings with 
non-equal but commensurable periods the technique we 
are about to describe uses the minimum common pe- 
riod for the two gratings). The medium in between is 
assumed to be homogeneous (it will become evident in 
the following that this condition can be easily relaxed). 
Both structures lie in the (x, y) plane of a orthonormal 
cartesian coordinate system. The z coordinate is called 
longitudinal because it is normal to both grating sur- 
faces. The choice of a rectangular basis is not manda- 
tory but dictated by symmetry, and it has been shown 
that sometimes the introduction of a nonrectangular co- 
ordinate system can improve the numerical convergence 

As usual for periodic system [49] [50l, it is conve- 
nient to introduce the transverse reciprocal space basis 
|k||, nm, A) associated with the reciprocal lattice de- 
scribed by the two dimensional vector 



_ 2nn ^ 2irm ^ 

\\,nm k|| + ~j— X + —j— y 

= a n x + /3 m y n,mGZ, (2) 
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where we have implicitly defined 

k|| = /c x x + k y y, 
«n = fcx + and p m = 



2nrn 



operator only connects wave vectors belonging to differ- 
ent Brillouin zones. Let us consider that the object "1" 
(3) is located to the left of the object "2" (i.e. if z\ < Z2 
since both gratings are parallel to the (x,y) plane). The 

(4) 

reflection operators take the form 



L x and L y are the periods of the grating along the x and 
y directions, respectively. The polarization is denoted 
by A, and is chosen to be one of the orthogonal polar- 
ization, s or p, where s-polarized fields have the electric 
field along the direction e s = Km nm x z/|K|| ?nm x z|, 
and p-polarized fields have the electric field along the di- 
rection e p = K|| 5nm x e s /|K|| 5nm x e s |. The transverse 
wavevector k|| is constrained by —tt/L x < k x < tt/L x 
and —Tr/Ly < k y < 7r/L y , which define the domain B 
known as the first Brillouin zone [49, 50 of the recipro- 
cal lattice (n = m = 0), the other Brillouin zones are 
defined by n, m ^ 0. The symbol ^ indicates forwards 
(—>) or backwards (<—) propagation. The correspond- 

(c) 

ing longitudinal wavevector ±/c^,nm is automatically ob- 
tained from the value of the frequency uj\ = i& and trans- 
verse wavevector, K||, together with the appropriate cav- 
ity material dispersion relation [51] 



eMi)% = K 2 



1 1 ,nm 



(5) 



Here Rek, 



(c) 



Im kifnm > and e c (uj) is the permit- 



tivity of the medium in the cavity formed by the two 
gratings. The two coordinate basis are indeed connected 
by a transformation which is equivalent to a three dimen- 
sional Fourier series and we have 



Si 



& 



(7) 



Denoting by O a generic block matrix, we then have 



(k||,nm ) A|Q|kf|,n'm',A'> = Q X ^, n , m ,5{^ -kf,). (8) 



Because of the symmetry, in this basis translation oper- 
ators behave in a very similar way 



X, 
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x< 



21 
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(9) 



where X are diagonal matrices with elements 



» _ 



| | ,nm 



•R-«V 



(10) 



-Ikiirn = je c (m?/C 2 + Kl nm - 



Using that "Trlog = logdet" and the Leibniz formula 
for the determinant of block matrices, Eq.Q takes the 
form [42] 



(r|k||, nm,^, A) e 



| | ,nm ' 



~R±ik 



(c) 



(6) 



where R = (x,y). The basis is also orthonormal. 

All operators in can now be represented as (infinite) 
matrices. Due to the periodicity and the symmetry, each 



Ho) = logdet 1 - * ^( a ) * ^2 • £(a) 



or more explicitly 



(ii) 



Ho) 



1=0 



d»ky 



log det 



c-A;A 7 



nm\n'm' 



£i(k|h^) 



A;A" 



nm:n" m' 



^2(k||,i6 



-cl(k, 



(c) 



(12) 



where the Einstein convention was used. Alternatively, 
one can also first diagonalize the matrix and then take 
the sum of the logarithm of each eigenvalue [52J. Eq. 
( 12 ) clearly reduces the problem of the calculation of the 



Casimir free energy to the determination of the matrix 
elements 7^(k||, ^i)nm-n'm" We will see in the next sec- 
tion that these quantities correspond to the (Rayleigh) 
reflection coefficients associated with the scattering of the 
electromagnetic field from an isolated grating. 



III. MODAL EXPANSION FOR 2D GRATINGS 

Several well developed theoretical and computational 
techniques exist for the evaluation of scattering from pe- 
riodic structures [22]. One class of scattering techniques 
makes use of the planewave expansion for the electromag- 
netic fields in the spirit of the derivation of the previous 
section. These techniques are typically referred to gener- 
ally as rigorous coupled wave approaches (RCWA) [25] . 
or specifically as Fourier modal techniques or planewave 
techniques. In this section, we will adapt thisl technique 
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FIG. 2: Schematic representation of a generic periodic struc- 
ture. Each unit cell can be thought of as the result of the 
superposition of different lamellar layers with the same pe- 
riod. This can be used to calculate scattering properties of 
arbitrary 3D scatterers provided the period is fixed. An inci- 
dent planewave in media is scattered into reflected diffration 
orders, and the transmitted diffraction orders in L+lth media 
are illustrated. 



to compute the Casimir force from a periodically pat- 
terned substrate. Each periodic structure will be thought 
of as the result of the superposition of 2D lamellar grat- 
ings (see Fig{2|. For each layer one can derive a direct 
eigenvalue problem for Maxwell's equations or the asso- 
ciated wave equation [48]. The corresponding complex 
eigenmodes in the 2D stratum will form a natural basis 
for the electromagnetic fields in the structure and also 
a representation of the scattered fields. The scattering 
properties of the total multilayer medium (lamellar or 
non lamellar) will be obtained by an iterative application 
of the scattering matrix following the theory of optical 
networks. To help the physical intuition in this section 
we work with real frequencies uj. However, at the end, 
the relevant quantities that enter in (12) are functions of 

We will 



the imaginary Matsubara frequencies uji 
discuss the analytical continuation of the modal approach 
to complex frequencies at the end of this section. 



A. The eigenvalue problem in the modulated 
region 

To begin with let us consider a single layer 2D peri- 
odic structure described by the complex dielectric per- 
mittivity e(x,y) and the complex magnetic permeability 
fi(x,y), both assumed to be scalars for simplicity. The 
permittivity and the permeability are periodic functions 
with period L x and L y of the transverse coordinates 
(x,y). Our eigenmode calculation starts by decompos- 
ing Maxwell's equations into transverse and longitudinal 
components and follows the approach of L. Li [20 , 48 , 53J. 
This splitting yields the so called waveguide equations 

- ikod z E t = V t [xz • V t x H t ] - feg/xz x H t , (13a) 



ikod z H t = -V t [Cz-VtxEt] 



•fcgezxEt, (13b) 

with x(x,y) = l/e(x,y), ((x,y) = l/fi(x,y), and 
fco = ^/ c - We assume that all fields depend on time 
as e~ lut . The longitudinal electric E z and magnetic H z 
fields are not independent and can be determined from 
the transverse components 



E z 



k e 



• V+ x H+ and H z 



z-V t xE £ . (14) 



Eqs. (13a) and (13b) further imply that the displace- 
ment field and magnetic induction are divergence free. 
We have not made any assu mptio ns re gardi ng the trans- 
verse mode structure in Eqs. (13a) and (13b) since in gen- 



eral the modes will not be classified as, e.g., transverse 
electric (TE) or transverse magnetic (TM) in a general 
complex metallic/dielectric composite structure. 

According to the Floquet-Bloch theorem inside each 
layer the field is a pseudoperiodic function and can be 
decomposed as follows 



f(r) = ^f„ m (z) e* K 



| | ,nm 



R 



(15) 



where we have combined the transverse components of 
the electric and magnetic fields into a single vector 



f(r) = 

where E. 



( E x (r) 
E y (r) 
H x (r) 

\ H v (t) 



fnm(^) — 



H V ' nm \ Z \ I • ( 16 ) 

- L1 x,nm / 



x(y),nm an d Ef x (y),nm are the Fourier coefficients 
of the transverse electric and magnetic fields, respec- 
tively. We similarly decompose the dielectric permittiv- 
ity, magnetic permeability and their inverses in Fourier 
series 



e(x,y) 


\ ^ i27rnx 1 L x -\-i2irmy / L y 

/ j tfim^ i 

nm 


(17) 


fi(x,y) 


\ ^ i27rnx / L x -\-i27rmy / L y 

— / j f^nm^- i 

nm 


(18) 




\ ^ i27rnx/ L x -\-i27rmy / L y 

— / Anm e ' 5 

nm 


(19) 




\ ^ /■ %2'Knx I 'L x -\-i2irmy / 'L y 

— / j snm& 

nm 


(20) 



5 



Collecting all the Fourier coefficients f nm in one single 
large vector F, it is possible to show that the waveguide 
equations (13) can be recast as a first order matrix dif- 



ferential equation [52j [53 



The matrix H has a block form and each bloch element 
is given by [48] 



■ik d z F(z)=H-F(z). 



(21) 






-cy 8 f C [nn ' ] 

>.[nn'] 



\ -R R /-l nn J , u£A nn J 






nn '] _ ^ 2 £ [ nn/ ] 

mm'] [mm'] 

R Ann'] 



a n >f3mX 



[nn ] 



[nn'] . ; 2 [nn 7 ] \ 
-anttn'Xjw] +^0M[W] \ 
/3 [nn 7 ] 



/ 



(22) 



where a n = fc^ + 2irn/L x and /3 m = k y + 2irm/L y . The 
symbol means that we have to take the shifted 

Fourier expansions of the complex permittivity, perme- 
ability, and their inverses (e.g. ej™^ = e n _ n /, m _ m /) that 
derive from the application of the Laurent rule [54] . The 
solution of the first order matrix differential equation ( |2l] ) 
has the form 



(23) 



where Y and 7 are respectively one of the eigenvectors 
and the corresponding eigenvalue that are solutions to 
the eigevalue problem [52] [55] 



n-Y. 



(24) 



Note that the value of the transverse momentum in the 
first Brillouin zone k|| is fixed in these equations, so that 
the set of eigenvectors Y u and corresponding eigenvalues 
7^ are functions of k|| . Given the fact that 1-L is non- 
hermitian, the eigenvalues are in general complex and 
the eigenvectors are not orthogonal [52j[55] to each other. 
However the eigenvectors are bi-orthogonal [52 , 55] to the 
eigenvectors of the corresponding adjoint equation 



AfegYt 



Y f -n ] . 



(25) 



From the theory of non-self adjoint differential equations 
one knows that the eigenvalues and eigenvectors of the 
two mutually adjoint equations can be ordered in such a 
way that \ v = 7^. Defining the scalar product 



we have that 



{Yl\Y v ,) = 8 vy . 
Note that the dimension of the matrix 1-L is even 
dim[H] =2D x 2D, 



(26) 



(27) 



(28) 



where D = 2(2N + 1)(2M + 1) and 2N + 1, 2M + 1 are 
the number of Fourier terms in the series expansion along 
the x and y directions, respectively. The characteristic 
equation for the eigenvalues is an equation of an even 
order with 2D solutions: if 7^ is a solution then 7* is 
also a solution, so that half of the eigenvalues and half of 
the eigenvectors can be obtained by simple conjugation. 
Notation-wise it is also possible to associate to the same 
value 7^ one eigenvector of eq.(21) and one eigenvector 
of the adjoint equation, i.e. 



^ Yjy, Yj , 



(29) 



and, in literature, Y v and Yj* are generally called right- 
handed and left-handed eigenvector respectively (the bi- 
orthogonality property can also be reinterpreted in terms 
of the previous definition). Furthermore, if we rearrange 
the vector F(z) to have the form 



F(z) 



_ fF E (z) 



(30) 



where Fe and Fh contains all electric field and magnetic 
field components, respectively.Then the matrix T-L takes 
the block form 



n 



H H 
He 



(31) 



The eigenvalue problem (24) can be reduced to separate 
equations [48] 



7^Y e =Hh-H e -Y e , 
A 2 Y h =He-%-Y h , 



(32a) 
(32b) 



so that we can deduce that the eigenvalues come in pairs 
and so do the eigenvectors 



and 



(33) 



representing forward and backward (<— ) propaga- 
tion, respectively. 
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FIG. 3: Schematic description of the scattering process by 
a generic object. The scattering matrix formalism connects 
the output field amplitudes as a function of the input field 
amplitudes for a general scattering object. 



These considerations ensure that the previous eigen- 
modes are natural basis for the electromagnetic field in 
the modulated region. Collecting all these results and 
definitions we have 



F(*) = ££"( k n) e " 



(k,i: 



B v . (34) 



Here we have explicitly indicated the dependency of the 
eigenvectors on k|| . The expression for the total field f (r) 
can be obtained by isolating the four dimensional vectors 
f nm from the previous expression and inserting them in 
Eq.(15). The only unknowns in Eq.p4|) are the coef- 
ficients A v and B V1 the field amplitudes (scalars), and 
they will be completely determined in the next subsec- 
tion once we impose the boundary conditions on the fields 
at the border of each modulated region. Note that the 
Rayleigh basis defined in Section [IT] for the cavity is a par- 
ticular (homogeneous medium / no modulation) case of 
the modal solution described above, and the eigenvalues 

(c) 

coincide with the longitudinal vector kz,hm, with mul- 
tiplicity two (one for each polarization A = s,p). This 
correspondence means that we can safely use the same 
formalism throughout the structure. 



B. Scattering and boundary conditions 

Let us now consider the problem of the scattering of an 
electromagnetic wave by a periodic structure in its com- 
pleteness, and assume that the network element of Figj3] 
is composed by several 2D periodic layers with thick- 
ness hi, electric permittivity 6i(x,y), and magnetic per- 
meability Hi(x, y) (i = 1 . . . L, see Figj2j. The first (1|0) 
interface is at the right of the network element at po- 
sition z = 0, while the last {L + 1\L) is at the left of 
the network at z = —h = Y^i^i- The first (eo,/io an d 
z > 0) and the last media (e^+ij/iL+i and z < —h) are 
assumed to be uniform layers with planewave type eigen- 
modes (Rayleigh basis). For each layer we can write the 



expression given in Eq.p4|) in the following form 

F (<) W=y (i) -4V-(B«) ' (35) 
The matrix y W contains all eigenvectors as columns 

y (i) = Ei, ■ • ■ , % { X D ) ■ (36) 



The block matrix 











(37) 



has the exponential functions with the eigenvalues 7^ 
along the diagonal of each of the two D x D block sub- 
matrices and describes the propagation through the layer 
of thickness hi . The field amplitudes for each layer A$ , 
By^ are collected in the column vectors and B^. 
That means that the vectors A^ ^ and B^ ^ with com- 
ponents A u ^ and B u ^ give the incident and reflected 
field amplitudes from the right of the network. Similarly 
A( L+1 ) and B( L+1 ) are the reflected and incident fields 
from the left of the network. The right and left reflection 
matrices are then defined as (see Fig. [3| 

B<°>=£.A<°>, A< L+1 ) = 5-B( i+1 ). (38) 

Similarly the left and right transmission operators can be 
defined as follows 



A (L+1) = • A (0) , B (0) = 7^ • B (L+1) . 



(39) 



Imposing the continuity of the tangential fields at the 
z = Zi = — Y^j=i hj interface means that F^ +1 ^(^) = 
F^(^) and the transfer operator that relates the field 
amplitudes in the ith and i + 1th layer is defined as 



= £ (i+i|<) . 



where 



t (i+i\i) = 



11 
21 



12 



(40) 



(41) 



Each single block of the transfer matrix has a dimension 
D x D. Defining rectangular D x 2D matrices 



1 



(42) 



the block elements can be expressed as overlaps of the 
mode eigenvectors 



t 



(i+lli) 

11 


= ^;(*+l)t 
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(i+lli) 

12 


= ^;(i+l)t 




(i+l|i) 
21 


= y(<+l)t 
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(43) 
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The field at the point z = Zi+i is related to that at z = Zi 
by 

F (i+1 \z i+1 ) = O^+^F^izi), (44) 

with 

e a+i\i) = t (i+m . t w (45) 

being the total transfer matrix from layer i to layer i + 1. 

Now we could construct the transfer matrix of the 
whole structure by iterating through the multilayer 
(lamellar or non-lamellar) structure, and then solve for 



the reflection operator. However, if numerically imple- 
mented, the transfer matrix method is known to suffer 
from instabilities when the layers are thick, due to the 
growing exponentials contained in the transfer matrix 
[56] . A remedy to these numerical instabilities is the 
S-matrix approach. The S-matrix is derived from the 
ordered T-matrix, where the ordering is in terms of for- 
ward and back propagating modes (see Fig{3|. The layer 
S-matrix can be defined by reordering the coefficients, 
and can be expressed in terms of the interface transfer 
matrices 



a+wi) a+iu) r n 1 Tj(i+i) = ^ • nu+i) ■ y^) 



(47) 



B W J ~ \ e^ (i)h > ) ' ^ ag +llt) J ' \ B^ 1 ) J = * '\ B^+D 

The interface S-matrix is 

[ Ai+m _(i+l|<) \ / Ai+Mi) _ Ai+l\i)rAi+l\ih-lAi+l\i) r.(i+l|i)i_i 

a ll a 12 1 — 6 H 6 12 L 6 22 J L 21 L 12 l L 22 J 

Ai+Mi) _(i+l|<) ~ r.(i+l\ih- lf (i+l\i) 

\ a 21 a 22 J \ ~l l 22 \ l 21 



The S-matrix stability is guaranteed in Eq. (47) due 
to the exponential decay (Inry* > 0) of the backward 
propagating mode. The layer S-matrix 



(<+i|i) (i+i|i) 

5 11 5 12 
*21 5 22 



(48) 



incorporates the propagation phase factor and the inter- 
face S-matrices. Again, the S-matrix of the multilayer 
object can be constructed using an iterative procedure 
with the following recursion relations 



V(*+1|0) _ Ai+l\i) fi Q (i+l|i)y(i|0)\ 1 v (i|0) 

^11 — s ll ^ — s 21 ^12 J ^11 j 

V (<+1|0) _ , A _ Q (i+l|i) v (i|0)\ _1 (i+l|i) y (i|0) 

^12 — 6 12 i" 6 H ^ b 21 n 12 J b 22 ^12 J 

y(*+l|0) _ v (i|0) v (i|0) A _ c (i+l|i) v (i|0)\ _1 (i+lK)y,(i|0) 

^21 — ^21 "T" ^22 ^ *21 ^12 y *21 ^11 > 

^22 +1 '°^ = -£ H1W Eg 0) (l - 4i +1 'Mf )" • (49) 



It is clear that the scattering matrix of the multilayer 
scatterer which connects the incident media (incident and 
reflected modes) to the exit media (transmitted modes) 
is given by £( L+1 I ). Unfortunately the definition of the 
scattering matrix in the grating literature [48] does not 
match the one used in the network theory (see Section 
|TT| ). However this last one can be obtained by a simple 
reordering of the block matrices 



From here all the relevant reflection (and transmission) 
operators for the calculation of the Casimir free energy 
can be obtained. 



C. Analytical properties of the modal approach 

In the previous section we derived the scattering ma- 
trix for a general 2D-periodic structure having in mind 
a real value for the frequency of the field. For the cal- 
culation of the Casimir free energy it is more convenient 
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to work with imaginary frequencies from the very begin- 
ning, since the quantities entering in Eq.(12) are func- 
tions of uj = z£. From an analysis of the previous section 
one can easily realize that the consequence of a (Wick) 
rotation in the complex plane (from the real frequency 
axis to the imaginary frequency axis) is automatically 
connected with the analytical properties of the involved 
functions. Maxwell waveguide equations (13) can be eas- 



ily written in terms of imaginary frequency. The dielec- 
tric permittivity and the magnetic permeability are an- 
alytic functions in upper part of the complex plane with 
the property [51 



e(z)=e*(-z*), t i(z)= f i*(-z*), 



(51) 



from which it immediately follows that e(z£) and /i(z£) are 
real. Since we are dealing with passive media they also 
are positive quantities for £ > 0. Clearly now the matrix 
1~L is real. By looking at the elements of 1-L one can easily 
show that its characteristic equation A2£>(z, ce, /?, 7 2 ) = 
inherits some of the properties of the permittivity and of 
the permeability, in particular it is true that 

A 2D (z, a, /3, 7 2 ) = A 2D (-z*, a, ft [ 7 2 ]*) (52) 

which means that the eigenvalues satisfy the property 
7^(z,a,/3) = z*, a, /?)]*. The choice of the defini- 

tion of the square root and its continuity in the complex 
plane finally implies that 



1 v {z,ol,P) = -7*(-**,a,/3), 



(53) 



from which it is possible to conclude that j u (u,a,f3) = 
— 7*(— co>, a, /?), and that 7zy (z£, a, /?) (£ G Re) is pure 
imaginary. Since they depend only upon the sign of the 
eigevalues, the concepts of forward and backward propa- 
gation can be easily generalized and we can still write an 
expansion like (34), working with imaginary frequency. 



The main difference is now that the exponentials are real 
functions. The formalism of the scattering matrix can 
be reapplied to derive the reflection operators at imag- 
inary frequencies. Naturally, in the numeric calculation 
the necessary truncation of the matrix T~L leads to some 
small deviations from the previous conclusions. For ex- 
ample the eigenvalues as a function of the imaginary fre- 
quency can acquire a small real part. 

In summary, to calculate the Casimir free energy be- 
tween periodic structures it is sufficient to write the 
waveguide equation for each Matsubara frequency and 
perform the corresponding modal expansion procedure. 
The resulting reflection matrices are then inserted in 
Eq.(12) and one can choose between the calculation of 



the determinant of the final expression, or its diagonal- 
ization followed by the calculation of the trace of the 
resulting diagonal matrix. Special attention is required 
when computing the zero Matsubara frequency contribu- 
tion to the free energy £ = in the scattering formalism. 
This will be outlined in the next section. 



D. Zero frequency modal solutions 

The computation of the Casimir free energy requires 
explicitly the I = (zero Matsubara frequency £o = 0) 
scattering matrices from the modal solution to Maxwell's 
equations. In the following we will consider non-magnetic 
media for simplicity, n{x,y) = 1. The zero frequency 
limit is intrinsically tied to dispersive models of the per- 
mittivity since in general the product uje(uj] x, y) appears 
explicitly in Maxwell's equations. Therefore, if the per- 
mittivity is finite at zero frequency then the electric and 
magnetic field are irrotational. However, if the permit- 
tivity has a simple pole at zero frequency the fields are 
coupled. 

The simplest and most useful model for the permittiv- 
ity of a homogeneous dielectric medium is the oscillator 
model or Drude-Lorentz model [51] . This model assumes 
harmonically bound charges to an ion core that makes 
up a neutralizing background. The resulting dielectric 
permittivity is 



e{u>) 



n 2 



(54) 



where Y is the damping coefficient, ujq is the oscillator 
frequency, and the plasma frequency is defined as 



O 2 

11 P l — 



47re 2 n 



m 



(55) 



This simple model of a dispersive media can be general- 
ized to multiple oscillators with different resonance fre- 
quencies, oscillator strengths and linewidths, and forms 
the basis for models of optical dispersion in dielectrics 
and metals. The Drude model for a metal is found by 
considering the limit of free electrons (i.e ujo — >• of Eq. 
([54]) ) and is given by 



c(oj) = 1 



Q 2 
lL pl 



iYuj' 



(56) 



The Drude model has two poles, one pole at zero fre- 
quency and the other pole in the lower half of the complex 
plane [28 . The connection with constituent relationships 
in Maxwell's equations is obtained by expanding Eq.(56), 



e(oj) = 1 



T(w + iT) Flo ' 



(57) 



and we can identify the first two terms as the Debye per- 
mittivity, and the last term with the Drude conductivity, 
a = Qpi/T. In the limit r — » 0, the Drude conductivity is 
infinite and we recover the London's superconductor [57] 
or also plasma model for the metallic media, which is 
singular at zero frequency. It is therefore important the 
order in which the limits are taken, and this can yield dif- 
fering results [3j[58]. In all our calculations, we consider 
the zero frequency limit first with finite conductivity at 
zero frequency. 
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Gasimir Force 



v J Distance (a) 



FIG. 4: AFM setup for measuring Casimir interactions be- 
tween a sphere and a ID lamellar grating. 



Th e zero fr equen cy limit to Maxwell's equations (Eq. 
( pa) and Eq.(|l3b|)) is 

c 



d z E t = V t 



z • V t x H t 



and 



9,Ht = VA-— zxEt, 



(58) 
(59) 



where we have used the Drude model of Eq.(|57|), and 
a = cr(x,y) is the spatially varying DC conductivity of 
the 2D periodic structure. These two equations replace 
the waveguide equations and are the basis for the zero 
frequency modal solutions. The corresponding transverse 
fields are the zero frequency limits of the transverse prop- 
agating solutions, and as such are different from purely 
static solutions. The zero frequency waveguide equations 
depend on the z component of the magnetic field. We 
therefore consider only transverse magnetic (TM) zero 
frequency solutions, H z = 0, since a non-zero value would 
imply a DC surface current is flowing, which is inconsis- 
tent with a non zero dissipation. Furthermore, for a pla- 
nar metallic mirror the zero frequency Fresnel reflection 
matrix for the transverse electric (TE) mode vanishes, 
and we are therefore restricted to TM solutions (H z — 0) 
that contribute to the zero Matsubara frequency in the 
expression for Casimir free energy. As a technical remark 
we point out that in the regions within the structure con- 
taining vacuum the conductivity a(x, y) is obviously zero, 
and the equations above are singular. In the numerics, 
we choose a vanishingly small but non-zero conductiv- 
ity for those regions. The final results turn out to be 
independent of this choice. Alternatively, we can con- 
sider small but finite frequencies approaching the zero 
frequency limit using the assumed material dispersion. 
Both approaches lead to the same converged zero fre- 
quency contribution for the grating structures examined. 



IV. RESULTS 

Most experiments measuring the Casimir force be- 
tween two objects use a sphere or a spherical lens as one 



of the intervening objects. This avoids parallelism issues 
which could affect the precision required by the measure- 
ments (see Fig. |4|. Although the scattering approach 
allows one to take into account spherical or even more 
complicate geometries [39] (47] 159H6T] , the modal method 
presented in this paper is not suited for non-periodic, 
spherical surfaces. However, the radius of curvature R of 
the sphere (or spherical lens) used in experiments is so 
much larger than the distance to the other plate that one 
can safely use the so-called proximity force approxima- 
tion (PFA) [62 , 63 . In our case PFA relates the Casimir 
force between a spherical surface and a grating to the free 
energy per unit of area between a plane and a grating 



F sph - g (a) = 2ttR T % 



pi- gi 



(60) 



and the force gradient in the sphere-grating configuration 
is then related to the force per unit area in the plane- 
grating geometry. The quantity on the r.h.s. can be 
calculated starting from the results presented in the pre- 
vious sections. In this section we benchmark our finite- 
temperature numerical code against a recent precise mea- 
surement of the Casimir force gradient between a metallic 
sphere and a deeply etched (« 1 micron) ID lamellar sil- 
icon grating [31]. We then move on to more complex 2D 
periodic structures. 

The ID lamellar grating in the experiment of Ref.[31 
was p-doped Si and the sphere was metallized with Au. 
In order to have a precise comparison between our numer- 
ics and the experimental data one should input into our 
code the actual optical properties of the samples used in 
the experiment. Unfortunately these were not measured 
in [31], so here we will use tabulated optical data for Au 
and p-doped Si, that have been compiled and studied by 
several authors [64H66] . We model the intrinsic Si per- 
mittivity by a Drude-Lorentz model, 



+ (e 



(61) 



with e = 11.87, = 1.035 , and u = 6.6 x 10 15 rad/s. 
The p-doped Si is modeled by adding to the intrinsic part 
a Drude background 



(62) 



with 



3.6151 x 10 14 rad/sec and 7 = 7.868 x 10 



13 



rad/sec. Similarly, the Au sphere is modeled by a Drude 
model 



* l p 



(63) 



with n p = 1.27524 x 10 16 rad/sec and T = 6.59631 x 10 13 
rad/sec. 

With the materials parameters at hand, we can be- 
gin comparing the experimental data and the computed 
exact scattering solution for the grating. Our code is 
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designed for 2D periodic systems, so in order to treat 
ID gratings we use a period along the translational in- 
variant direction (say y) equal to that along the non- 
invariant direction x, i.e. L x = L y , so that the resulting 
structure is 2D periodic. In all our calculations we will 
truncate all the Fourier series to N = M = 5. This 
implies that the dimension of the reflection matrices is 
2(2N + l) 2 x 2(2N + l) 2 = 242 x 242. The k||-space 
integration is performed using a 16pt Gauss-Legendre 
quadrature. We set the temperature to T = 300K, and 
use the first 36 Matsubara frequencies for all the studied 
range of distances between the sphere and the grating. 

We consider the two ID lamellar gratings used in 
[31]. The Casimir force in the plane-grating configura- 
tion is evaluated by numerical differentiation of the plane- 
grating free energy. The resulting force gradient in the 
sphere-gradient geometry calculated with our numerical 
code is shown in Figs. [5] and [6] where we also show the 
experimental data with their errors, kindly provided to 
us by H.B. Chan. Given the uncertainty in the optical 
parameters used in our numerics compared those of the 
actual samples, the experiment-theory agreement seems 
to be very reasonable. A calculation of the reduced chi 
square gives Xred ~ 2.9 for sample A and Xred ~ 8. 8 for 
sample B. One can also test the deviations of the exact 
numerical results from the prediction of the proximity 
force approximation (PFA). In this approximation, the 
force between the plane and the deeply etched grating 
is computed by multiplying the force between two plane 
(i.e., the usual Lifshitz force for the plane-plane geome- 
try) by the filling factor of the grating (we neglect the 
contribution of the bottom part of the grating to the 
PFA result since its depth is sufficiently large). In the 
insets of Figs. [5] and [6] we plot the ratio between the 
exact and PFA results for the Casimir force gradients in 
the sphere-grating configuration. There are three sets 
of data represented in those insets: a) the solid line is 
the ratio of our exact numerics and the theoretical PFA 
prediction based on Lifshitz theory, both using the pa- 
rameters above for the optical data of the samples; b) 
the hollow circles with their error bars are the ratio of 
the experimental data of [31 for the force gradient in 
sphere-grating geometry with respect to the "experimen- 
tal" PFA. This last data set was obtained from a separate 
measurement of the force gradient in the sphere-plane 
geometry and subsequently multiplied by the grating's 
filling factor [31]. Note that since the plane and the grat- 
ing in [31] were fabricated following identical procedures, 
one can expect them to have the same optical properties; 
c) the squares represent the ratio of our exact numerics 
and the "experimental" PFA. Given the uncertainty in 
the optical parameters used in the numerics, the most 
unbiased comparison is to compare case b) against case 
c). In that case, one is effectively comparing numerators 
normalized by the same denominator. Again, the theory- 
experiment comparison is reasonably good in view of the 
experimental and theoretical uncertainties. 

The numerical errors in the computation of the Casimir 
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FIG. 5: Casimir force gradient in the sphere-grating geom- 
etry for sample A in [31]: exact numerics (solid) and ex- 
perimental data (dots with error bars). The inset shows 
the ratio of the Casimir force gradient divided by the PFA 
prediction: exact numerics/theoretical PFA (solid), exact 
numerics/ "experimental" PFA (squares), and experimental 
data/ "experimental" PFA (hollow circles with error bars). 
The geometrical parameters of grating A are: period=1000 
nm, depth=980 nm, and filling factor=0.510. The radius of 
the sphere is 50/rni, and the temperature is set to T = 300K. 
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FIG. 6: Casimir force gradient in the sphere-grating geom- 
etry for sample B in [31]: exact numerics (solid) and ex- 
perimental data (dots with error bars). The inset shows 
the ratio of the Casimir force gradient divided by the PFA 
prediction: exact numerics/theoretical PFA (solid), exact 
numerics/ "experimental" PFA (squares), and experimental 
data/ "experimental" PFA (hollow circles with error bars). 
The geometrical parameters of grating B are: period=400 
nm, depth=1070 nm, and filling factor=0.478. The radius of 
the sphere is 50/xm, and the temperature is set to T = 300K. 



force originate from three sources: 1) the truncation of 
the Matsubara frequency summation, 2) the truncation 
of the discrete spatial frequency spectrum resulting in 
finite dimensional reflection matrices, and 3) numerical 
integration over the continuous transverse wavevector in 
the first Brillioun zone. The truncation of the Matsub- 
ara summation is determined by the temperature and 
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FIG. 7: Computed Casimir force gradient for ID and 2D Si 
gratings. Together with the ID grating of sample B (dot 
dashed curve) two simple but representative 2D extensions 
have been considered. The first example is an array of p- 
doped Si pillars, and the second is the complementary struc- 
ture, i.e. a free standing membrane etched with an array of 
square holes. The period of the structures is 400 nm and 
the filling factor of the pillar and of the membrane are 1/4 
and 3/4, respectively. The etch depth for all the structures is 
1070 nm. For all cases, a Au sphere with radius R = 50/mi 
has been used. It has also been assumed that the tempera- 
ture is T = 300K. We find that the force gradient scales as 
the filling factor. The deviation from the PFA is shown in 
the inset and scales this time with inverse of the filling factor, 
having the strongest deviation in the case of the pillars. 



the minimum distance required in the force displacement 
curve. This minimum distance sets the maximum Mat- 
subara frequency for fixed temperature, such that for a 
minimum spacing of 100 nm at 300K, we find that 36 
Matsubara frequencies give free energy convergence of 
better than 10 -4 . This frequency cutoff is used in all 
numerical results presented here. 

The truncation of the discrete spatial frequency 
(diffraction orders) is necessitated by the need to deal 
with finite matrices. Our modal approach improves with 
increasing number of diffraction orders but so does the 
computational cost. By increasing the number of diffrac- 
tion orders, we find that the short distance behaviour for 
the grating is improved. We have used N = 5 for all 
calculations, which was observed to be enough to ensure 
the convergence at approximately 1% accuracy. It is im- 
portant to note that matrices are dense and scale as 47V 4 , 
which limits practically the spatial frequency cutoff. 

The final source of numerical error is the numerical in- 
tegration over the continuous transverse wavevector in 
the first Brillioun zone. It is desirable to have mini- 
mal /c-space sampling to reduce the assembly and com- 
putation of the reflection matrices for each Matsubara 
frequency. The integration of the continuous wavevec- 
tor is performed numerically using a an n-point Gauss 
Legendre quadrature. In order to estimate the integra- 
tion error, we have computed the the force gradient for 5 



different n-point (9,16,25,36,49) Gaussian quadratures in 
the 1st Brillioun zone. The relative error, the standard 
deviation divided by the mean, for the experimental dis- 
placement range is less than 3% and quite small at small 
displacement. This is quite reasonable since the force 
gradient at large separations is quite small so small stan- 
dard deviations give large relative error. From our three 
estimates of error, our estimated total error in computing 
the Casimir force is expected to be less than 3% over the 
entire force gradient displacement curve. 

The detailed analysis of the ID lamellar grating prob- 
lem allowed us to validate the scattering approach by 
comparing to high precision experimental data. Here we 
want to extend the method to 2D periodic structures 
which provide more tailorable properties in which we can 
manipulate the Casimir force. We will consider two sim- 
ple but representative extensions of the ID grating. The 
first example is an array of p-doped Si pillars. Simi- 
lar geometries have been used to obtain a negative in- 
dex of refraction [57] in the optical range [68 , 69 . They 
have also been exploited to investigate the phenomenon 
of quantum reflection [70] [71] of an atom over the purely 
attractive Casimir-Polder potential generated by the pe- 
riodic structure. The second case is the complementary 
structure, i.e. a free standing membrane etched with an 
array of square holes. This structure is very similar to the 
one used to measure the extraordinary light transmission 
through sub- wavelength apertures [72H74] . 

In Fig (7] we show the prediction for Casimir force gra- 
dient for the sample B ID grating (dash-dotted line), the 
pillars (dotted line), and the membrane (dashed line). 
The period of the structures is in all cases 400 nm and 
the eteched depth is 1070nm, and the filling factor of the 
pillars and of the membrane are 1/4 and 3/4, respec- 
tively. We find that the Casimir force gradient scales as 
the filling factor, i.e. the force is less attractive for the 
pillars than for the membrane, with the case of grating 
in between these two cases. The inset shows the compar- 
ison with the relative PFA approximation (filling factor 
times the Lifshitz force between plane and sphere). If the 
PFA were valid, the ratio (inset) would be unity for all 
separations. The deviation from the additive approxima- 
tion scales with the inverse of the filling factor, having 
the strongest deviation in the case of the pillars. 



V. DISCUSSION AND CONCLUSION 

We have developed a finite-temperature modal ap- 
proach to compute Casimir interactions between 2D pe- 
riodic structures. We have compared our computational 
approach to high precision published experimental data 
for ID lamellar gratings. This benchmark validated our 
modal approach and led to good agreement between the- 
ory and experiment, as confirmed by a reduced chi square 
values of 2.9 and 8.8 for the two samples used in [31]. 

In order to demonstrate the flexibility of our approach, 
we also calculated the Casimir force between the first 
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simple extensions of a ID grating, namely an array of 
square pillars and an array of square holes. In both cases, 
as already known for the ID grating, we observe deviation 
of the Casimir force gradient from the value obtained 
from the Proximity Force Approximation. This deviation 
scales with the filling factor and it is more accentuated 
in the case of the pillars. 

We plan in the near future to extend these results to 
complex metallic structures, such as 3D metamaterials, 
which were recognized as possible candidates to engineer 
the Casimir force between two vacuum separated objects. 
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